Mon. Not. R. Astron. Soc. 000, [TJjse] (2008) Printed 30 April 2009 (MN WT^ style file v2.2) 



Scaling for the intensity of radiation in spherical and 
aspherical planetary nebulae 

L. Zaninetti 

Dipartimento di Fisica Generate, 

Via Pietro Giuria 1 

10125 Torino, Italy 

email: zaninetti@ph.unito.it 

to be inserted 

ABSTRACT 

The image of planetary nebulae is made by three different physical processes. 
The first process is the expansion of the shell that can be modeled by the 
canonical laws of motion in the spherical case and by the momentum conser- 
vation when gradients of density are present in the interstellar medium. The 
second process is the diffusion of particles that radiate from the advancing 
layer. The 3D diffusion from a sphere as well as the ID diffusion with drift 
are analyzed. The third process is the composition of the image through an 
integral operation along the line of sight. The developed framework is applied 
to A39 , to the Ring nebula and to the etched hourglass nebula MyCn 18. 

Key words: ISM: jets and outflows , ISM: kinematics and dynamics , ISM: 
lines and bands , planetary nebulae: individual 



1 INTRODUCTION 



The planetary nebula , in the following PN , rarely presents a circular shape generally 
thought to be the projection of a sphere on the sky. In order to explain the properties of PN, 



Kwok et al. (1978) proposed the interacting stellar wind (ISW) theory. Later on Sabbadin 



et al. (1984) proposed the two wind model and the two phase model. More often various 



types of shapes such as elliptical , bipolar or cigar are present, see Balick (1987); Schwarz 



et al. (1992); Manchado et al. (1996); Guerrero et al. (2004); Soker fc Hadar (2002); Soker 



(2002). The bipolar PNs , for example , are explained by the interaction of the winds which 
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originate from the central star , see Icke (1988); Frank et al. (1995); Langer et al. (1999); 



Gonzalez et al. (2004). Another class of models explains some basic structures in PNs through 



hydrodynamical models, see Kahn & West (1985); Mellema et al. (1991) or through self- 



organized magnetohydrodynamic (MHD) plasma configurations with radial flow, see Tsui 



(2008). 



An attempt to make a catalog of line profiles using various shapes observed in real PNs 



was done by Morisset & Stasinska (2008). This ONLINE atlas , available at 



http://132.248.1.102/Atlas_profiles/img/, is composed of 26 photo-ionization models corre- 



sponding to 5 geometries, 3 angular density laws and 2 cavity sizes, four velocity fields for 
a total of 104 PNs, each of which can be observed from 3 different directions. 



Matsumoto et al. (2006) suggest that a planetary nebula is formed and evolves by the 



interaction of a fast wind from a central star with a slow wind from its progenitor , an 
Asymptotic Giant Branch (AGB) star. It seems therefore reasonable to assume that the PN 
evolves in a previously ejected medium ( AGB) phase in which density is considerably higher 
than the interstellar medium (ISM) . We can , for example , consider a PN resulting from a 
5 Mq Main Sequence (MS) star . The central core will be a White Dwarf (WD) less than 1 
Mq and the ionized nebula is generally less than 1 Mq . We therefore have ^ 3 Mq of gas 
around the PN which come from the AGB. The number density that characterizes the PN 
is 

9.66M5^ particles 



n 



R3 



cm-' 



(1) 



where M^^ q is the number of solar masses in the volume occupied by the nebula and Rpc 
the radius of the nebula in pc. 



By inserting M0=O.6O5, see for example Figure 2 in Perinotto et al. (2004) , and Rpc=l 



in the previous formula we obtain n ^ 6.28 



particles 



. This can be considered an averaged value 



and it should be noted that the various hydrodynamical models give densities ,p , that scale 
with the distance from the center R as R~" , with 2.5 < a < 3.5 , see Villaver et al. (2002), 



Perinotto et al. (2004), Schonberner et al. (2005), Schonberner et al. (2005), Schonberner 



et al. (2007), and Steffen et al. (2008). 



The already cited models concerning the PNs leave a series of questions unanswered or 
partially answered: 

• Which are the laws of motion that regulate the expansion of PN ? 

• Is it possible to build up a diffusive model in the thick advancing layer ? 
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• Is it possible to deduce some analytical formulas for the intensity profiles ? 

In order to answer these questions Section |2] describes three observed morphologies of 
PNs, Section[3]analyzes three different laws of motion that model the spherical and aspherical 
expansion, Section |4] reviews old and new formulae on diffusion and Section |5] contains 
detailed information on how to build an image of a PN. 



2 THREE MORPHOLOGICAL TYPES OF PNS 

This section presents the astronomical data of a nearly spherical PN known as A39, a weakly 
asymmetric shell , the Ring nebula , and a bipolar PN which is the etched hourglass nebula 
MyCn 18 . 



2.1 A circular spherical PN 

The PN A39 is extremely round and therefore can be considered an example of spherical 



symmetry, see for example Figure 1 in Jacoby et al. (2001 ) . In A39 the radius of the shell , 

Rshell is 

R,hcii = 2.42 X 10^^677021 cm = 0.78 pc , (2) 
where B77 is the angular radius in units of 77" and D21 the distance in units of 2.1 kpc , 



see 



Jacoby et al. (2001) . The expansion velocity has a range [32 ^ 37 ^] according to 



Hippelein & Weinberger (1990) and the age of the free expansion is 23000 yr, see Jacoby 



et al. (2001). The angular thickness of the shell is 
(5rsheii = 3.17 lO^^GioDsicm = 0.103 pc , 



(3) 



where Gio is the thickness in units of 10.1" and the height above the galactic plane is 1.42 



kpc , see Jacoby et al. (2001). The radial distribution of the intensity in [QUI] image of A39 
after subtracting the contribution of the central star is well described by a spherical shell 
with a 10" rim thickness, see Figure p!) and [Jacoby et al.| ( |200l| . 

The caption of Figure [l] also reports the of the fit computed according to formula 
In presence of real data a merit function , is introduced as 



X 



^Yf-'^? , (4) 

where N is the number of the data , yi the theoretical ith point , yi,obs the ith observed point 
and (Ji the error for the ith observed point here computed as ^. 
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Arcsec from Nebula Center 



Figure 1. Cut of the relative intensity of A39 crossing the center in the east- west direction (dotted Une with some error bar) 
and the rim model (full line) fully described in Jacoby et al. (2001) , = 0.505 when 5819 point are considered. 



2.2 The asymmetric PN 



The Ring nebula , also known as M57 or NGC6720 , presents an elliptical shape characterized 
by a semi-major axis of 42", a semi-minor axis of 29.4" and ellipticity of 0.7, see Table I in 



Hiriart (2004). The distance of the Ring nebula is not very well known ; according to Harris 



et al. (1997) the distance is 705 pc . In physical units the two radii are 

Rsheii,minor = O.I629.4D705 pc scmi — miuor radius 
Rsheii.major = O.I4642D705 pc scmi - major radius , 



(5) 



where 629.4 is the angular minor radius in units of 29.4", G42 is the angular major radius 
in units of 42" and D705 the distance in units of 705 pc. The radial velocity structure in the 
Ring Nebula was derived from observations of the H2 (molecular Hydrogen) v = 1- S(l) 
emission line at 2.122 /xm obtained by using a cooled Fabry- Perot etalon and a near- infrared 



imaging detector , see Hiriart (|2004|) . The velocity structure of the Ring Nebula covers the 
range [-30.3 ^ 48.8 



km l 
s J 



2.3 The case of MyCn 18 

MyCn 18 is a PN at a distance of 2.4 kpc and clearly shows an hourglass-shaped nebula, see 



Corradi & Schwarz 


(1993 


); 


Sahai et al. 


(1999 



we can fix the equatorial radius in 2.80 x 10^^ cm , or 0.09 pc , and the radius at 60° from 
the equatorial plane 3.16 x 10^^ cm or 0.102 pc . The determination of the observed field of 



velocity of MyCn 18 varies from an overall value of 10 — as suggested by the expansion of 



[OIII] , see Sahai et al. (1999) , to a theoretical model by Dayal et al. (2000) in which the 
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velocity is 9.6 ^ when the latitude is ° (equatorial plane) to 40.9 ^ when the latitude is 
60 °. 



3 LAW OF MOTION 

This Section presents two solutions for the law of motion that describe asymmetric ex- 
pansion. The momentum conservation is then applied in cases where the density of the 
interstellar medium is not constant but regulated by exponential behavior. 



3.1 Spherical Symmetry - Sedov solution 

The momentum conservation is applied to a conical section of radius R with a solid angle 



-(AMR) = AF^ 



A f2, in polar coordinates, see McCray & Layzer (1987) 
d 
dt 

where 

e 

10 



(6) 



AM 



p{R,e,<f))dV, (7) 

is the mass of swept-up interstellar medium in the solid angle A Q, p the density of the 
medium , P the interior pressure and the driving force: 

A F = FR'^An . (8) 



After some algebra the Sedov solution is obtained, see Sedov (1959); McCray & Layzer 



( [1987| ) 

R(t) = 



25 Et 



2\ 1/5 



V, 4 vr p y 

where E is the energy injected in the process and t the time. 



(9) 



Another slightly different solution is formula (7.56) in Dyson & Williams (1997) 



R(t) 



'25 Et 



2\ 1/5 



(10) 



^3 TT p y 

where the difference is due to the adopted approximations. 

Our astrophysical units are: time (t4), which is expressed in 10^ yr units; E42, the energy 
in 10'^^ erg; and uq the number density expressed in particles cm~^ (density p = ngm, where 
m=1.4mH). With these units equation ^ becomes 



R(t) ^ 0.198 



E^2 t4 

. no 



2\ 1/5 



pc 



(11) 



6 L. Zaninetti 

The expansion velocity is 

which expressed in astrophysical units is 
V(t) ^ 7.746 



E, 



km 



(12) 



(13) 



By inserting M0=O.6O5 and Rpc=l in formula Q we obtain n ^ 6.282S£|^ . This value is 



higher than the value of number density of the ISM at the plane of the galaxy, n ^ 1 



particles 



. Equations (11) and (13) represent a system of two equations in two unknowns : t4 and E42 



. By inserting for example R = 0.78 pc in equation (11 ) we find 
t4 = 77.15 ^ 



E 



■42 



and inserting V = 35 km s ^ in equation (13) we obtain 
0.3954/^^ = 35 . 



(14) 



(15) 



The previous equation is solved for E42 = 7833.4 that according to equation (14) means 



t4=.87173. These two parameters allows a rough evaluation of the mechanical luminosity 
L = Y that turns out to be L ^ 2.847 lO^^ergs s~^. This value should be bigger than the 
observed luminosities in the various bands. As an example the X-ray luminosity of PNs , 
Lx, in the wavelength band 5-28 A has a range [10'^°'^ ^ lO'^^'^ergs s^^] , see Table 3 in 



Steffen et al. (2008) 



Due to the fact that is difficult to compute the volume in an asymmetric expansion the 
Sedov solution is adopted only in this paragraph. 



3.2 Spherical Symmetry - Momentum Conservation 

The thin layer approximation assumes that all the swept-up gas accumulates infinitely in a 
thin shell just after the shock front. The conservation of the radial momentum requires that 

^7rRVR = Mo , (16) 

where R and R are the radius and the velocity of the advancing shock , p the density of the 
ambient medium , Mq the momentum evaluated at t = to , Ro the initial radius and Rq the 



initial velocity , see Dyson & Williams (1997); Padmanabhan (2001). The law of motion is 

1 

Ro _ " 



R = Ro(l + 4^(t-to) 
V ^ 



(17) 



Scaling for the intensity of radiation in spherical and aspherical planetary nebulae 7 
and the velocity 

_ 3 

R = Rofl + 4|^(t-to)l ' . (18) 



R, 







From equation (17) we can extract Rq and insert it in equation (18) 
''-4(t-t„) RJ V^-RT) ■ 

The astrophysical units are: t4 and to,4 which are t and to expressed in 10"^ yr units, Rpc and 
Ro,pc which are R and Rq expressed in pc, Rkms and Ro,kms which are R and Rq expressed in 
Therefore the previous formula becomes 

_ 3 

Ri.m. = 24.49 \ .^V^^' ^+ ^'tT^^"" ) ' • (20) 

On introducing Ro,pc = 0.1 , Rpc = 0.78 , Rkms = 34. , the approximated age of A39 is 
found to be t4 — to,4 = 50 and Ro,kms = 181.2. 

3.3 Asymmetry - Momentum Conservation 

Given the Cartesian coordinate system (x, y, z) , the plane z = will be called equatorial 
plane and in polar coordinates z = Rsin(6'), where 6 is the polar angle and R the distance 
from the origin . The presence of a non homogeneous medium in which the expansion takes 
place can be modeled assuming an exponential behavior for the number of particles of the 
type 

z R X sin(6') . ^, 

n(z)=noexp-- =noexp , (21) 

where R is the radius of the shell , uq is the number of particles at R = Rq and h the scale. 
The 3D expansion will be characterized by the following properties 

• Dependence of the momentary radius of the shell on the polar angle 6 that has a range 
[-90° ^ +90°]. 

• Independence of the momentary radius of the shell from (p , the azimuthal angle in the 
x-y plane, that has a range [0° ^ 360°] . 

The mass swept, M, along the solid angle A Q, between and R is 

M(R) = ^mHnoI„,(R) + ^vrR^nomn , (22) 
where 

U(R)= / r^exp -^dr , (23) 
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where Rq is the initial radius and nin the mass of the hydrogen . The integral is 

h f 2 h2 + 2 Roh sin (6) + Rq^ (sin (9)?) e" 
Im(R) = ^ — ^ 



Rq sin(9) 
h 



Rsin(fl) 
h 



h (2 h^ + 2 Rh sin (9) + R^ (sin {9)f) e 

3 • (2^) 



(sin(0)) 

The conservation of the momentum gives 

M(R)R = M(Ro)Ro , (25) 

where R is the velocity at R and Rq the initial velocity at R = Rq. 

In this differential equation of the first order in R the variable can be separated and the 
integration term by term gives 

M(r)dr = M(Ro)Ro x (t - to) , (26) 

where t is the time and to the time at Ro- The resulting non linear equation JF^l expressed 
in astrophysical units is 

^NL = -6e V hpc - hpce V (sin(6')) Ro,pc 

_ Ro,pc Bin(9) _ Rp.pc sin(fl) 

-e/ip/e V sin (9) Ro,pc- 3 hp,^e v {sin {9)y Ro,pc^ 

Rpc sin(fl) B.pc Bin(9) 

- Ro,pc (sin {9)r + 6 e v /i^/ + 4 e v hp^^Rpc sin (9) 

_ Rpc si.i(9) „ „ „ 

+e V hpc^Rpc^sin{9)y 

+2e V /ip//? 

pc sin (^) + 2e /ip^ (sin(^)) Ro,pc 

+e V /ipc-Rpc (sin(6')) /?o,pc 
+ (sin {9)f Ro,pc^Rpc - 0.01 (sin (9)^ i?o,pc'Ro,kms (t4 - to,4) = , (27) 

where t4 and to,4 are t and to expressed in 10^ yr units, Rpc and Ro,pc are R and Rq expressed 
in pc, Rkms and Ro.kms are R and Ro expressed in 6* is expressed in radians and hp^ is 
the the scale , h , expressed in pc. It is not possible to find Rpc analytically and a numerical 
method should be implemented. In our case in order to find the root of J?-nl, the FORTRAN 



SUBROUTINE ZRIDDR from Press et al. (1992) has been used 



The unknown parameter t4 — to,4 can be found from different runs of the code once Ro,pc 
is fixed as ~ 1/10 of the observed equatorial radius , Ro,kms is 200 or less and hpc ~ 2 x Ro,pc. 

From a practical point of view, e , the percentage of reliability of our code can also be 
introduced, 

e = (1 - l(Rpc,°bs-Rpc.num)| ^ _ ^ ^28) 

Rpc,obs 
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Table 1 . Data of the simulation of the Ring nebula 

Initial expansion velocity ,Ro,kms 200 

Age(t4-to,4) ' 0.12 

Initial radius Ro,pc 0.035 

scaUng hpc 2 X Ro,pc 




Figure 2. Continuous three-dimensional surface of the Ring nebula : the three Eulerian angles characterizing the point of 
view are $=180 °, 0=90 ° and *=-30 °. Physical parameters as in Table [T] 

where Rpc,obs is the radius as given by the astronomical observations in parsec , and Rpc,num 
the radius obtained from our simulation in parsec. 

In order to test the simulation over different angles, an observational percentage of reli- 
ability ,eobs, is introduced which uses both the size and the shape, 

eobs = 100(1 - ^j|Rpc.°bs-Rpc.num|j ^^ ^29) 
Z^j -tlpCjObs j 

where the index j varies from 1 to the number of available observations. 



3.3.1 Simulation of the Ring nebula 

A typical set of parameters that allows us to simulate the Ring nebula is reported in Table [T} 
The complex 3D behavior of the advancing Ring nebula is reported in Figure |2] and 
Figure [3] reports the asymmetric expansion in a section crossing the center. In order to 
better visualize the asymmetries Figure |4] and Figure |5] report the radius and the velocity 
as a function of the position angle 6. The combined effect of spatial asymmetry and field of 
velocity are reported in Figure |6j 
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-0.1 -aoB 



X fpc] 



Figure 3. Section of the Ring nebula on the x-z plane. The horizontal and vertical axis are in pc. Physical parameters as in 
Table [l] 




9 [degree] 

Figure 4. Radius in pc of the Ring nebula as a function of the position angle in degrees. Physical parameters as in Table [T] 



The efficiency of our code in reproducing the observed radii as given by formula (28 
and the efficiency when the age is five time greater are reported in Table [2] 

An analogous formula allows us to compute the efficiency in the computation of the 
maximum velocity , see Table |3j 



3.3.2 Simulation of MyCn 18 

A typical set of parameters that allows us to simulate MyCn 18 is reported in Table |4| 
The bipolar behavior of the advancing MyCn 18 is reported in Figure [7] and Figure [8] 
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Figure 5. Velocity in ^ of the Ring nebula as a function of the position angle in degrees. Physical parameters as in Table^ 




Figure 6. Map of the expansion velocity in relative to the simulation of the Ring nebula when 300000 random points are 
selected on the surface. Physical parameters as in Table [T] 



Table 2. Reliability of the radii of the Ring nebula. 





Rup{pc) polar direction 


Rcq{pc) equatorial plane 


Robs 


0.14 


0.1 


Rnum(our code) 


0.125 


0.102 




89 


97 


e (%) for a time 5 times greater 


27 


41 
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Table 3. Reliability of the velocity of the Ring nebula 



V( — ) maximum velocity 



Vobs 48.79 

Vnum 39.43 

e(%) 80.81 

e(%) (%) for a time 5 times greater 35.67 



Table 4. Data of the simulation of MyCn 18 



Initial expansion velocity ,Ro,kms [km s ^] 200 
Age (t4-to,4) [10-* yr] 0.2 
Initial radius Ro,pc [pc] 0.001 
scaling h [pc] 1.0 X Rq 



reports the expansion in a section crossing the center. It is interesting to point out the 



similarities between our Figure of MyCn 18 and Figure 1 in Morisset & Stasinska (2008) 



which define the parameters a and h of the Atlas of synthetic line profiles. In order to better 
visualize the two lobes Figure [9] reports the radius as a function of the position angle 6. 



The combined effect of spatial asymmetry and field of velocity are reported in Figure 10 



The efficiency of our code in reproducing the spatial shape over 12 directions of MyCn 



18 as given by formula (29 ) is reported in Table |5j This Table also reports the efficiency in 
simulating the shape of the velocity. 



Figure 11 reports our results as well those of Table 1 in Dayal et al. (2000). 




Figure 7. Continuous three-dimensional surface of MyCn 18 : the three Eulerian angles characterizing the point of view are 
$=130 °, 0=40 ° and 1'=5 °. Physical parameters as in Table [I] 



Scaling for the intensity of radiation in spherical and aspherical planetary nebulae 13 




Figure 8. Section of MyCn 18 on the x-z plane. Physical parameters as in Tablejl] 



Table 5. Reliability of the spatial and velocity shape of MyCn 18. 



radius velocity 
«obs(%) 90.66 57.68 



4 DIFFUSION 

The mathematical diffusion allows us to follow the number density of particles from high 
values (injection) to low values (absorption). We recall that the number density is expressed 
in p-^rticies symbol C is used in the mathematical diffusion and the symbol n in an 

unit volume 
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Figure 9. Radius in pc of MyCn 18 as a function of latitude from 0° to 60° ( dotted line) when the physical parameters are 
those of Table|4] The points with error bar (1/10 of the value) represent the data of Table 1 in Dayal et al. 2000. 




Figure 10. Map of the expansion velocity in relative to the simulation of MyCn 18 when 300000 random points are 
selected on the surface. Physical parameters as in Table |4] 
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6 [degree] 



Figure 11. Velocity in of MyCn 18 as a function of the latitude in degrees when the physical parameters are those of 
Table|4] dotted line. The points with error bar (1/10 of the value) represent the data of Table 1 in Dayal et al. 2000. 

astrophysical context. The density p is obtained by multiplying n by the mass of hydrogen 



mn , and by a multiplicative factor , f, which varies from 1.27 in Kim et al. (2000) to 1.4 



m 



McCray & Layzer (1987) 



p = fmnn 



(30) 



The physical process that allows the particles to diffuse is hidden in the mathematical dif- 
fusion. In our case the physical process can be the random walk with a time step equal to 
the Larmor gyroradius. In the Monte Carlo diffusion the step-length of the random walk 
is generally taken as a fraction of the side of the considered box. Both mathematical diffu- 
sion and Monte Carlo diffusion use the concept of absorbing-boundary which is the spatial 
coordinate where the diffusion path terminates. 

In the following, 3D mathematical diffusion from a sphere and ID mathematical as well 
Monte Carlo diffusion in presence of drift are considered. 

4.1 3D diffusion from a spherical source 

Once the number density , C, and the diffusion coefficient ,D, are introduced , Fick' s first 
equation changes expression on the basis of the adopted environment , see for example 



equation (2.5) in Berg (1993). In three dimensions it is 

dC 



dt 



DV^C 



(31) 



where t is the time and is the Laplacian differential operator. 
In presence of the steady state condition: 



DVX = 



(32) 
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The number density rises from at r=a to a maximum value Cm at r=b and then falls 



again to at r=c . The solution to equation (32) is 
r 

where A and B are determined by the boundary conditions 



Cab(r) = Cm 1 



and 



Cbc(r) = a 



- 1 



- 1 



a ^ r ^ b 



b < r < c 



Berg 


(1993 


) or in 


Crank 


(1979 



(33) 



(34) 



(35) 



4.2 ID diffusion with drift, mathematical diffusion 

In one dimension and in the presence of a drift velocity ,u, along the radial direction the 



diffusion is governed by Pick's second equation , see equation (4.5) in Berg (1993) , 



dC _ ^d^C ^dC 
~dt ~ 



(36) 



where vi can take two directions. The number density rises from at r=a to a maximum 



value Cm at r=b and then falls again to at r=c . The general solution to equation (36) in 
presence of a steady state is 



C(r) = A + Bei 



(37) 



We now assume that u and r do not have the same direction and therefore u is negative ; 
the solution is 



C(r) = A + Be-E-- , 

and now the velocity u is a scalar. 
The boundary-conditions give 



(38) 



Ca,b,drift(r) — C^ 

and 

Cb,c,drift(r) = Cj] 



e d" — e D 

i^ii Tih a ^ r ^ b downstream side , 

e — e D 



b ^ r ^ c upstream side . 



(39) 



(40) 



e — e D ^ 

A typical plot of the number density for different values of the diffusion coefficient is reported 
in Figure [121 
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pc from Nebula Center 



Figure 12. Number density of A39 as a function of the distance in pc from the injection when u = 1 , Cm = 1, a = 69.6 arcsec, 
b = 87 arcsec , c = 89.6 arcsec and D = 2 (full line ), D = 7 (dashed ), D = 12 (dot-dash-dot-dash) and D = 17 (dotted ). The 
conversion from arcsec to pc is done assuming a distance of 2100 pc for A39. 



4.3 ID diffusion with drift, random walk 

Given a ID segment of length side we can implement the random walk with step-length A 
by introducing the numerical parameter NDIM = ^ . We now report the adopted rules 
when the injection is in the middle of the grid : 



(i) The first of the NPART particles is chosen. 

(ii) The random walk of a particle starts in the middle of the grid. The probabilities of 
having one step are pi in the negative direction (downstream) ,pi = | — x |, and p2 in the 
positive direction (upstream) , p2 = | + /.f x |, where /i is a parameter that characterizes 
the asymmetry (0 ^ // ^ 1). 

(iii) When the particle reaches one of the two absorbing points , the motion starts another 
time from (ii) with a different diffusing pattern. 

(iv) The number of visits is recorded on Ai , a one-dimensional grid. 

(v) The random walk terminates when all the NPART particles are processed. 

(vi) For the sake of normalization the one-dimensional visitation or number density grid 
M is divided by NPART. 

There is a systematic change of the average particle position along the x-direction: 

(dx)=A^A , (41) 

for each time step. If the time step is dt = where Vtr is the transport velocity, the 
asymmetry , that characterizes the random walk is 
u 

Vtr 



(42) 
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Figure 13. Number density in A39 of the ID asymmetric random walk (full line), NDIM=401 ,NPART=200 ,side = 40 arcsec 
, A = 0.1 arcsec and fj, =- 0.013. For astrophysical purposes /i is negative. The theoretical number density as represented by 
formulas ( |39[ l and ( |40[ l is reported when u = 1 , Cm = 1, a = 60 arcsec, b = 80 arcsec , c = 100 arcsec and D = 3.84 (dotted 
line ). The conversion from arcsec to pc is done assuming a distance of 2100 pc for A39. 



Figure 13 reports M.{x), the number of visits generated by the Monte Carlo simulation as 



well as the mathematical solution represented by formulas (39) and (40). 



The solutions of the mathematical diffusion equations (39) and (40) can be rewritten at 
the light of the random walk and are 



Ca,b,Mc(l') 

and 

Cb,c,Mc(r) 



_ 

e A — e 



2ii 



2m, 



c 



_ 2m 

e >■ 

m _2m„ 

e >- 



_2Mb 

e A " 



_2m 

e A 

„-2iib 



a ^ r < b downstream side 



b ^ r ^ c upstream side 



(43) 



(44) 



5 THE IMAGE OF THE PN 

The image of a PN can be easily modeled once an analytical or numerical law for the intensity 
of emission as a function of the radial distance from the center is given. Simple analytical 
results for the radial intensity can be deduced in the rim model when the length of the layer 
and the number density are constants and in the spherical model when the number density 
is constant. 

The integration of the solutions to the mathematical diffusion along the line of sight 
allows us to deduce analytical formulas in the spherical case. The complexity of the intensity 
in the aspherical case can be attached only from a numerical point of view. 
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5.1 Radiative transfer equation 

The transfer equation in the presence of emission only , see for example Rybicki & Lightman 



(1985 


) or 


Hjellming 


(1988) 



ds 



(45) 



where li, is the specific intensity , s is the line of sight , the emission coefficient, k^, a 
mass absorption coefficient, ( the mass density at position s and the index i> denotes the 



interested frequency of emission. The solution to equation (45) is 



(46) 



where Tj, is the optical depth at frequency u 

dr, = k,Cds . (47) 

We now continue analyzing the case of an optically thin layer in which r^, is very small ( or 
kj, very small ) and the density ( is substituted with our number density C(s) of particles. 
Two cases are taken into account : the emissivity is proportional to the number density and 
the emissivity is proportional to the square of the number density . In the linear case 

j,C = KC(s) , (48) 

where K is a constant function. This can be the case of synchrotron radiation from an 



ensemble of particles , see formula (1.175 ) in Lang (1999) . This non thermal radiation 
continuum emission was detected in a PN associated with a very long-period OH/IR variable 



star (V1018 Sco), see [Cohen et aL] ( |2006| ) . 
In the quadratic case 

j,C = K2C(s)2 , 

where K2 is a constant function. This is true for 



(49) 



Free-free radiation from a thermal plasma, see formula (1.219) in Lang (1999) . This 



radiation process was adopted by Gonzalez et al. (2006) in the httle Homunculus. 



Thermal bremsstrahlung and recombination radiation , see formula (1.237) in Lang 



(1999) . This radiation process was adopted in PNs by Blagrave "eFal] ( [2006| ); [Schwarz fc 



Monteiro (2006); Gruenwald & Aleman (2007). 



I.(s 



The intensity is now 

K / C(s/)ds/ optically thin layer linear case 

Jsn 



(50) 
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Figure 14. The two circles (section of spheres) which include the region with constant density are represented through a full 
line. The observer is situated along the x direction, and three lines of sight are indicated. 



or 

Iy(s) = K2 / C(s/)^ds/ optically thin layer quadratic case . (51) 

J So 

In the Monte Carlo experiments the number density is memorized on the grid Ai and the 
intensity is 

Hhi) = ^^^^ ^ -^ihh^) optically thin layer linear case , (52) 
k 

or 

/(i, j) = ^ A s X A^(i, j, k)^ optically thin layer quadratic case , (53) 
k 

where As is the spatial interval between the various values and the sum is performed over the 
interval of existence of the index k. The theoretical intensity is then obtained by integrating 
the intensity at a given frequency over the solid angle of the source. 



5.2 3D Constant Number density in a rim model 

We assume that the number density C is constant and in particular rises from at r = a 
to a maximum value Cm , remains constant up to r = b and then falls again to 0. This 



geometrical description is reported in Figure 14 The length of sight , when the observer is 
situated at the infinity of the x-axis , is the locus parallel to the x-axis which crosses the 
position y in a Cartesian x — y plane and terminates at the external circle of radius b. The 
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Figure 15. Cut of the mathematical intensity / of the rim model ( equation | |55| l) crossing the center (full line ) of A39 and 
real data (dotted line with some error bar ) . The number of data is 801 and for this model = 1-487 against = 0.862 of 
the rim model fully described in Jacoby et al. (2001). 



Table 6. Simulation of A39 with the rim model 



symbol meaning value 

a radius of the internal sphere 72.5" 

b radius of the external sphere 90.18" 

R-shell observed radius of the shell 77" 

'5rshell,t theoretical thickness of the shell 17.6" 

i^rshcll observed thickness of the shell 10.1" 

ratio of observed intensities (1.88 — 2.62) 

iccntcr ^ ^ 

T?'"— n> ratio of theoretical intensities 3.03 



locus length is 



loa = 2x (7b2-y2-^a2-y2) ;0^y<a 



U = 2x(^b2-y2) ;a^y<b . (54) 

When the number density Cm is constant between two spheres of radius a and b the intensity 
of radiation is 



Ioa = CmX2x (\/b2-y2_y^a2-y2) ;0<y 



< a 



U = CmX2x(Jb2-y2) ;a^y<b . (55) 



The comparison of observed data of A39 and the theoretical intensity is reported in Figure 15 
when data from Table |6] are used. 

The ratio between the theoretical intensity at the maximum , (y = b) , and at the 
minimum , (y = 0) , is given by 

l(y - b) 
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Figure 16. The spherical source inserted in the great box is represented through a dashed line, and the two absorbing 
boundaries with a full line. The observer is situated along the x direction, and three lines of sight are indicated. Adapted from 
Figure 3.1 by Berg (1993) . 

5.3 3D Constant Number density in a spherical model 

We assume that the number density C is constant in a sphere of radius a and then falls to 
0. 

The length of sight , when the observer is situated at the infinity of the x-axis , is the 
locus parallel to the x-axis which crosses the position y in a Cartesian x — y plane and 
terminates at the external circle of radius a. The locus length is 

U = 2x (y'a^-yS) ;0^y<a . (57) 

When the number density Cm is constant in the sphere of radius a the intensity of radiation 
is 

loa = Cm X 2 X (^a^-yS) ; ^ y < a . (58) 



5.4 3D diffusion from a sphere 



Figure 16 shows a spherical shell source of radius b between a spherical absorber of radius 
a and a spherical absorber of radius c. 

The number density rises from at r=a to a maximum value Cm at r=b and then falls 
again to at r=c . 



The numbers density to be used are formulas (34) and (35) once r = a/x^ + y^ is im 



posed ; these two numbers density are inserted in formula ( 49 ) which represents the transfer 
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equation with a quadratic dependence on the number density. An analogous case was solved 
Zaninetti (2007) by adopting a linear dependence on the number density . The geometry 



of the phenomena fixes three different zones (0 — a, a — b, b — c) in the variable y, ; the first 
piece , I^(y) , is 



.2 ,,2 



CmV , ^ 2 . .Va" -y 



>(b2-2ba + a2)(c2-2cb + b2)L ^cxc..™^ ^ 



2 a^ arctan( )cb — 2Ja? — y^ycb 



—2 ay \yi[J a? — y^ + a)b^ + 2 a^ arctan( )cb 



+2 yy^b2 - y2cb + 2 aln(y^b2 - y^ + b)yb^ + 2 aln(yb2 - y^ + b)yc 
-2 cy ln(y^b2 -y2 + b)b2 - 2 c^ arctan(^^^5!Ez!)ba - 2 y^b^^^ba - 



2 cy \n{^h'^ - y^ + h)a? + 2 c ln(yc2 - y^ + c)yb 



+2 c ln(Y/c2 - y2 + c)ya^ + 2^c^ - y^yba + 2 c^ arctan(^^ ^)ba 



y 



2ayln(\/a2 — y^ + a)c^ + ya^ — y^yc^ 



—4 c ln( Y c^ — y2 + c)yba + 4 ay ln( y a^ — y^ + a)cb 



— y^ya^ + a^ arctan( )b^ — yyb^ — y^c^ — c^ arctan( )s? 

— a/c^ — y2yb^ + y-y/b^ — y^a^ + c^ arctan( ^ )b^ 



—a? arctan( ^ '^^ )b^ + \J — y^yb^ 
— c^ arctanf ^ )b^ + a? arctanf ^ ^ )c^] (59) 

y y 

^ y < a . 

The second piece , I^^y) , is 



I"(y)= / 2CLdx+ r 2CLdx 



2 



■y(b2 - 2 ba + a2)(c2 - 2 cb + b^) [^V^^ " ^^^^ + arctan(— — )b 

— c^ arctan( ^ ^ )b^ 

b^ — y^a^ + c^ arctan( )a^ + c^ arctan( )b^ + J — y^ya^ 



-y\ - , , 

y y 



+ - y^yb^ + 2 a ln(y)yb2 + 2 a ln(y)yc2 + 2 cy \nUh^ - y^ + b)b2 
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—2 arctan( )cb — 2 y yb^ — y^cb — 2 aln(Yb^ — y^ + b)yb^ 

■ Vc^ - y^ ^ 



— 2 aln(Yb^ — + b)yc —2c arctan( )ba — 2 y — y^yba 

—2 c ln(y — y^ + c)ya^ — 2 c ln( y — y^ + c)yb^ + 2 arctan( )ba 

2 ' 



+2 yi/b2 - y2ba + 2 cy ln(i/b2 - y^ + b)a^ - 4 a ln(y)ycb + 4 c ln(^c2 - y2 + c)yba] (60) 

a ^ y < b . 

The third piece , l"^(y) , is 



F(y)=/^^"'^2CLdx 

2,2 . /V^C^ 

■ /, a , ON / a , , ox 1 7 1/ - , ^ + c arctaiil 

y(b2-2ba + a2)(c2-2cb + b2)^-' V v ^ 

~ y^yb^ + arctan( ^ )a^ + y^ c2 — y2a^ — \j a? — y^yc^ 



1^2/^ 2 ^ , /p2 _ ,^2 

2 ,, \,7, fyv/c2-y2b^ + c^arctan(^ 



+a^ arctan( )b^ + A/b2 — y2yc^ — -i/b2 — y2ya^ — c^ arctan( )b^ 

y y 

— a^ arctanf )b^ — a^ arctanf )c^ + 2 cy ln(wb2 — y2 + b)a^ 

y y ^ 

Va^ - y^ . 



+2 ay ln(y a2 — y2 + a)c +2a arctan( )cb — 2 aln(y b2 — y2 + b)yc 

Vb^^, 



-2 aln(^b2 - y2 + b)yb^ - 2 ^b2 - y2ycb - 2 a^ arctan(^^ ^)cb 



+2 cy ln(^b2 - y2 + b)b^ -2c ln(^c2 - y2 + c)yb^ -2c ln( ^c2 - y2 + c)ya^ 



2y\/c2 — y2ba — 2c arctan( )ba + 2ayln(A/a2 — y2 + a)b 



2 v/b^^ 



+2 wa2 — y2ycb + 2 d\P' — y2yba + 2 c arctan( )ba 



+4 c ln(^c2 - y2 + c)yba - 4 ay ln(^a2 - y2 + a)cb] (61) 

b ^ y < c . 



The profile of / made by the three pieces ( 60), ( 61 ) and (62), can be cahbrated on the 
real data of A39 and an acceptable match is realized adopting the parameters reported in 
Table [3 

The theoretical intensity can therefore be plotted as a function of the distance from the 



center , see Figure [17] or as an image , see Figure [18| 

The effect of the insertion of a threshold intensity , Itr, given by the observational tech- 
niques , is now analyzed. The threshold intensity can be parametrized to Imax, the maximum 
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Table 7. Simulation of A39 with 3D diffusion 

symbol meaning value 

a radius of the internal absorbing sphere 65.96" 

b radius of the shock 80" 

c radius of the external absorbing sphere 103.5" 

Rsholl observed radius of the shell 77" 
■^rshcll observed thickness of the shell 10.1" 

ratio of observed intensities (1.88 — 2.62) 

iccntcr ^ ^ 

T?'"— n> ratio of theoretical intensities 2.84 




Arcsec from Fjebula Center 



Figure 17. Cut of the mathematical intensity / (formulas ( |60| l, ( |61| and ( |62[ |) , crossing the center (full line ) of A39 and 
real data (dotted line with some error bar). The number of data is 801 and for this model = 19.03 against = 12.60 of 
the rim model fully described in Jacoby et al. (2001). 
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Figure 18. Contour map of / particularized to simulate A39. 
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-100 



Arcsec from Nebula center 



Figure 19. The same as Figure [Ts] but with Itr = Imax/2 



value of intensity characterizing the ring: a typical image with a hole is visible in Figure 19 

when Itr = Imax/2. 

The position of the minimum of / is at y = and the position of the maximum is situated 
at y = b. 

The ratio between the theoretical intensity at maximum , Imax at y = b , and at the 
minimum (y = 0) is given by 

(62) 



Imax Numerator 
I(y = 0) Denominator 
where 



Numerator = (b^ — 2 ba + a^ 

X (2 cb ln(b) -2c ln(Vc2 - b2 + c)b + hVc^ - + arctan( ^)) 

b 

and 



(63) 



Denominator = 

2 b(a^c - c^a - 2 bca ln(a) + 2 bca ln(c) - ba^ + bc^ - b^c + b^a + b^a ln(a) 
-c^aln(b) + b^cln(b) - b^aln(b) - b^cln(c) + a^cln(b) + c^aln(a) - a^cln(c)) . (64) 

The ratio rim(maximum) /center (minimum) of the observed intensities as well as the 
theoretical one are reported in Table [7j 



Up to now we have not described the fainter halo of A39 which according to Jacoby et al. 



(2001) extends 15" beyond the rim. The halo intensity can be modeled by introducing two 



different processes of diffusion characterized by different geometrical situations . The first 
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Figure 20. Cut of the mathematical intensity / that characterizes the two-phase diffusion ( full line ) and real data of A39 
(dotted line with some error bar). The number of data is 801 and for this model = 2.29 against = 12.606 of the rim 
model fully described in Jacoby et al. (2001). 
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Figure 21. Contour map of the decimal logarithm of / of the two-phase diffusion relative to A39. 



is represented by / made by the three pieces ( 60), ( 61) and (62), the second one is the 
intensity between a larger sphere ( r = 2 x c) and smaller sphere ( r = b) with constant 



density , see formula (55 



I = I(Cm,i5 ai, bi, ci, y) + Ioa(Cm,2, a2, ba, y) 



(65) 



where the numbers 1 and 2 stand for first process and second process . The second process 
with constant density will be characterized by a larger volume of the considered bigger sphere 
and smaller number density , i.e. Cni,2 Cm,i • A typical result of this two phase process 



is plotted in Figure |20] and the image reported in Figure [21] ; the adopted parameters are 
reported in Table [8] 
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Table 8. Simulation of A39 , halo comprised 

symbol meaning value 

ai radius of the internal absorbing sphere 66.3" 

bi radius of the shock 80" 

ci radius of the external absorbing sphere 103.5" 

Cm,i maximum numberdensity main diffusion 1 

a2 internal radius of the halo process 80" 

b2 external radius of the halo process 207" 

Cm. 2 numberdensity halo 0.045 

ratio of simulated intensities 2.85 




-0^ 0^ 

pc from Nebula Center 



Figure 22. Cut of the numerical intensity / crossing the center (dotted line ) of A39 when the drift is considered and real 
data (full line). The parameters are u = 1 , a = 69.6 arcsec, b = 87 arcsec , c = 89 arcsec and D = 6. In this case j^^J^qj =2.92 
. The number of data is 801 and for this model = 20.96 against = 10.36 of the rim model fully described in Jacoby et 
al. (2001). The conversion from arcsec to pc is done assuming a distance of 2100 pc for A39. 



5.5 3D diffusion from a sphere with drift 

The influence of advection on diffusion can be explored assuming that in 3D the number 
density scales in the radial direction as does the ID solution with drift represented by 



formulas (39) and (40) . This is an approximation due to the absence of Pick's second 



equation in 3D. Also here the geometry of the phenomena fixes three different zones (0 — 



a, a — b, b — c) in the variable y , see Figure 16 , and the intensity along the line of sight can 



be found by imposing r = a/x^ + y^ . In this case, the integral operation of the square of the 



number density which gives the intensity can be performed only numerically , see Figure 22 



5.6 3D complex morphologies 

The numerical approach to the intensity map can be implemented when the ellipsoid that 
characterizes the expansion surface of the PN has a constant thickness expressed , for ex- 
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ample , as rmm/f where is the minimum radius of the elhpsoid and f an integer. We 
remember that f = 12 has a physical basis in the symmetrical case , see IMcCray & Layzer 



(1987). The numerical algorithm that allows us to build the image is now outlined 



• A memory grid A^(i,j,k) that contains NDIM^ pixels is considered 

• The points of the thick ellipsoid are memorized on Ai 

• Each point of Ai has spatial coordinates x, y, z which can be represented by the following 
1x3 matrix ,A, 



A 



X 



(66) 



The point of view of the observer is characterized by the Eulerian angles (<I',G,\1/) and 



therefore by a total rotation 3x3 matrix , E , see Goldstein et al. (2002). The matrix point 



is now represented by the following 1x3 matrix , B, 



B = E- A 



(67) 



• The map in intensity is obtained by summing the points of the rotated images along a 
direction , for example along z , ( sum over the range of one index, for example k ). 



Figure 23 reports the rotated image of the Ring nebula and Figure 24 reports two cuts along 
the polar and equatorial directions. 



Figure 25 reports the comparison between a theoretical and observed east-west cut in 



that cross the center of the nebula, see Figure 1 in Garnett & Dinerstein (2001). 



A comparison can be made with the color composite image of Doppler-shifted H2 emission 



as represented in Figure 2 in Hiriart (2004) 



In order to explain some of the morphologies which characterize the PN's we first map 



MyCn 18 with the polar axis in the vertical direction , see map in intensity in Figure 26 The 



vertical and horizontal cut in intensity are reported in Figure 28 The point of view of the 
observer as modeled by the Euler angles increases the complexity of the shapes : Figure 27 



reports the after rotation image and Figure 29 the vertical and horizontal rotated cut. The 
after rotation image contains the double ring and an enhancement in intensity of the central 
region which characterize MyCn 18. 
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Figure 23. Map of the theoretical intensity of the Ring nebula. Physical parameters as in Table [l] and f=12 . The three 
Eulerian angles characterizing the point of view are <I>=180 °, ©=90 ° and ^=-30 °. 




Figure 24. Two cut of the mathematical intensity / crossing the center of the Ring nebula: equatorial cut (full lino) and polar 
cut (dotted line) . 

This central enhancement can be considered one of the various morphologies that the 
PNs present and is similar to model BLi — F in Figure 3 of the Atlas of synthetic line profiles 



by Morisset & Stasinska (2008). 



6 CONCLUSIONS 

Law of motion The law of motion in the case of a symmetric motion can be modeled 
by the Sedov Solution or the radial momentum conservation. These two models allow to 
determine the approximate age of A39 which is 8710 yr for the Sedov solution and 50000 yr 
for the radial momentum conservation. In presence of gradients as given , for example , by an 
exponential behavior, the solution is deduced through the radial momentum conservation. 
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Figure 25. Cut of the mathematical intensity / of the Ring Nebula crossing the center (full line ) and real data of (dotted 
line with some error bar ) . The number of data is 250 and for this model = 15.53 . The real data are extracted by the 
author from Figure 1 of Garnett and Dinerstein 2001. 
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Figure 26. Map of the theoretical intensity of MyCn 18 . Physical parameters as in Table |4] and f=12 . The three Eulerian 
angles characterizing the point of view are $=180 ° , ©=90 ° and \E'=0 °. 

The comparison with the astronomical data is now more comphcated and the single and 
multiple efficiency in the radius determination have been introduced. When , for example , 
MyCn 18 is considered , the multiple efficiency over 18 directions is 90.66% when the age of 
2000 yr is adopted. 
Diffusion 

The number density in a thick layer surrounding the ellipsoid of expansion can be con- 
sidered constant or variable from a maximum value to a minimum value with the growing 
or diminishing radius in respect to the expansion position. In the case of a variable number 



density the framework of the mathematical diffusion has been adopted, see formulas (34) 
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Figure 27. Map of the theoretical intensity of the rotated MyCn 18 . Physical parameters as in Table |4] and f=12 . The three 
Eulerian angles characterizing the point of view are <&=130 °, 0=40 ° and ^=5 °. 




Figure 28. Two cut of the mathematical intensity / crossing the center of MyCn 18 : equatorial cut (full line) and polar cut 
(dotted line) . Parameters as in Figure [26] 




Figure 29. Two cut of the mathematical intensity / crossing the center of the rotated MyCn 18 nebula; equatorial cut (full 
line) and polar cut (dotted line) . Parameters as in Figure [27] 
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and (35). The case of diffusion with drift has been analytically solved , see formulas (39) 



and (40), and the theoretical formulas have been compared with values generated by Monte 
Carlo simulations. 
Images 

The intensity of the image of a PN when the intensity of emission is proportional to the 
square of the number density can be computed through 

• an analytical evaluation of lines of sight when the number density is constant between 



two spheres or in one sphere , see formula (55) and formula (58). 

• analytical evaluation of integrals when the number density is variable , see formulas 



60), ( 61) and (62), when the motion is symmetric. In this framework it is also possible to 



build a two-phase diffusion model that allows us to reproduce the faint extended halo, see 



formula (65). 

• a numerical evaluation of integrals when the number density is variable , the drift is 



present and the motion is symmetric, see Section 5.5 



• a numerical evaluation of lines of sight when the motion is asymmetric, see Section 5^ 

In the case of A39 the of comparison between the theoretical and observed cut in 
intensity can be evaluated , see Table [9j From a careful evaluation of Table [9] it is possible 
to conclude that the models here considered produce which are slightly bigger than the 
rim model of Jacoby et al. (2001). When , conversely , a diffuse halo is considered the is 
smaller ; this makes the study of the interaction between PN and the surrounding halo an 
interesting field of research. 

Next step 

Here we have explored the conservation of the radial momentum in a medium with 
an exponential behavior of the type p oc exp _ ^^^^"(^) -which is symmetric in respect to 
the plane z = . The next target can be the analysis of the conservation of the radial 
momentum in a spherical symmetry of the type p oc R~" , see Section [ij In this case the 
spatially asymmetric motion can be obtained by considering the conservation of the radial 
momentum in a medium with density of the type p oc R~" x exp 
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Table 9. Data of the simulation of the Ring nebula 



model 

rim with fixed thickness 1.487 

difTusion 19.03 

difTusion with drift 20.96 

diffusion + halo 2.29 



X^ of reference, Jacoby et al. 2001 



0.862 
12.60 
10.36 
12.606 
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